Non-linear dynamics of a continuous spring–block model of earthquake faults

The continuous one-dimensional Burridge–Knopoff model is generalized by introducing plastic creep in addition to rigid sliding. The resulting equations, for an order parameter (sliding rate) and a control parameter (driving force), exhibit a velocity-strengthening and a velocity-softening instability. In the former regime, considered to be the analog of self-organized criticality in continuum systems, anomalous diffusion is described by a non-linear diffusion equation. The latter regime, characteristic of deterministic chaos, is described by a time-dependent Ginzburg–Landau equation. PACS numbers: 05.45.+b, 46.30.Pa, 82.40.Bj, 91.30.Px Typeset using REVTEX 1 (September 18, 2000) The dynamics of non-trivial spring-block models of earthquake faults has regained attention since Bak and Tang [1] argued that the earth crust may be considered to be in a self-organized critical state. They showed that if stick–slip dynamics of a slowly driven system is computed by a cellular-automaton algorithm the system evolves into a statistical steady state with power law correlations in time and in space. These correlations compare favorably to experimentally observed power law correlations for earthquakes (for example, to the Gutenberg–Richter law for the frequency distribution of energy release, or to the Omori law for the distribution of aftershocks [2]). The statistical steady state that exemplifies scaling properties similar to those at a critical point but does not require careful tuning of an external variable has been identified as self-organized criticality (SOC) [3–5]. A complementary approach to obtaining power law scaling has been based on continuous, deterministically chaotic models. In particular, the properties of continuous velocitysoftening models, whose simplest form is known as the Burridge–Knopoff model [6], have been studied extensively; see, for example Refs. [7,8]. These investigations showed that power law scaling is reproduced even in the absence of an explicit stochastic element in the continuous differential equation: randomness is introduced only via the initial conditions. Even though both approaches reproduce power-law scaling, the connection between continuum and cellular automata models has not been fully elucidated, nor has the origin of SOC been unequivocally identified. SOC in continuum systems has been attributed to singular diffusion [9] and to the presence of conservation laws that lead to a non-linear diffusion equation [10]. Recently, an attempt has been made to reconcile discrete (cellular automata) and continuum SOC versions by Gil and Sornette [11]. They proposed a Landau–Ginzburg theory of self-organized criticality based on the coupled dynamics of a control and an order parameter. SOC was identified with an uphill diffusion of the control parameter close to a critical point (spinodal point). However, they did not explicitly construct a physical model that exhibits the proposed feedback mechanism, and they found that external noise is essential in obtaining scaling. We present a generalization of the continuous, one-dimensional uniform Burridge– 2 (September 18, 2000) Knopoff (BK) model and we analyze its stability properties. The generalization is motivated by physical arguments. Since the model is general enough it is suggested that it may become a vehicle for the identification of similarities and differences between SOC and deterministic chaos. We show that depending on driving conditions the system either evolves into a self-organized critical state or its time evolution becomes deterministically chaotic. The standard BK model of earthquake faults describes the dynamics of displacements u of a slowly driven spring–block chain of masses m in the presence of a non-linear dynamic friction Φ. The masses (of characteristic extension ξ) are longitudinally coupled by coil springs of stiffness kl and they are transversally coupled to the driving interface by leaf springs of stiffness kt. The essential feature of the model is the choice of the friction force that generates a velocity-softening instability. The generalization we propose is to account for fault creep via the introduction of an internal variable. Specifically, under the imposed driving the fault responds by a combination of rigid translation us (sliding) and plastic displacement up (irreversible deformation by creep of some boundary layer of the fault). The corresponding equations of motion become (primes denote spatial derivatives and dots time derivatives) müs = F − Φ(u̇s + u̇p) , (1a) Ḟ = klξ u̇ s − kt(u̇s − v̄)− F − Fy τp . (1b) The first equation describes balance of forces (Newton’s law) and the second the time evolution of the driving shear force F . The external loading rate is ktv̄ where v̄ is the tectonic drift velocity. The first term on the RHS of Eq. 1b arises from longitudinal compression or tension, the second from elastic shear, and the last gives the force relaxation due to plastic deformation under the assumption that the plastic displacement rate depends linearly on the driving force, u̇p = F − Fy ktτp . (2) 3 (September 18, 2000) The critical force corresponding to plastic yielding is denoted by Fy and τp is the characteristic time of plastic relaxation of the shear forces. In what follows we neglect Fy, without loss of generality. In analogy to critical phenomena the two coupled differential equations may be interpreted as a system of equations that describe the dynamics of a control parameter F (or the plastic displacement rate u̇p) and an order parameter u̇s (the sliding rate). This compares with the proposed feedback mechanism of Gil and Sornette [11]. Moreover, such a description of the dynamics of the fault is reminiscent of previous analyses of laboratory friction data [12,13] where a set of two coupled differential equations for the friction stress, a parameter that characterizes the evolving state of the surface, and a constitutive equation have been proposed. For the constitutive relation of Ref. [14] and in the absence of inhomogeneities (u̇ s = 0) these equation may be re-written in a form similar to ours: the equation for the driving force 1b becomes identical, whereas the equation corresponding to Eq. 1a contains additional non-linearities whose physical origin is difficult to justify [15]. A schematic diagram of the friction force we will be considering is shown in Fig. 1: for small u̇ = u̇s + u̇p slow, stable creep is described whereas at higher values velocity-softening is obtained. Note that plastic deformation introduces memory-dependent friction which is related to the aging of the fault by plastic accommodation of its interface, while previous analyses of the BK model have usually considered only the velocity-softening part of the friction force. The importance of using a more realistic friction force that allows for viscous creep has already been noted by Carlson and Langer [16]. Time and space are rendered dimensionsless by defining t̃ = t/τp and x̃ = (kt/kl) x/ξ, while the dimensionless sliding rate is defined by e = u̇sτp/lp, the dimensionless shear force by f = F/(ktlp), and the dimensionless friction force by φ = Φ/(ktlp). By this scaling we have introduced a characteristic length scale lp = τpv̄. In terms of the scaled fields Eqs. 1 become ǫė = f − φ[v̄(e+ f)] , (3a) 4 (September 18, 2000) ḟ = e − e− f + 1 . (3b) Besides the set of constants defining the friction force φ Eqs. 3 depend on two parameters, the drift velocity v̄ and a new parameter ǫ = 1 τp ( m kt )1/2 , (4) which is the ratio of the natural frequency of transversal oscillations of individual blocks to the plastic relaxation time. It is tempting to identify the ǫ → 0 (τp → ∞) limit as the standard, homogeneous BK model. However, inspection of the scaling relations and Eq. 2 shows that the BK limit corresponds to τp → ∞ with u̇p → 0, namely creep is suppressed while the friction force remains constant. The natural limit of Eqs. 3 is ǫ → 0 keeping the plastic displacement rate u̇p fixed since the plastic relaxation time τp and the characteristic length scale lp go to infinity at the same rate (for fixed v̄). The physical interpretation of this limit is that an increase of the plastic relaxation time is associated with an increase of the driving force F since force relaxation by creep becomes less efficient. Consequently, the limit that gives the standard BK model is not physical reasonable in our framework since the shear force and the plastic relaxation time are interrelated (cf. Eq. 2) [15]. The limit ǫ → 0 corresponds to an adiabatic elimination of the fast variable e (the order parameter u̇s). This yields the evolution equation for the slow variable f (control parameter F ) ḟ = [D(f)f ] − 1 v̄ ζ(f) + 1 , (5a) D = 1 v̄ ∂fζ(f)− 1 , (5b) where we have defined ζ to be the inverse function of φ, ζ(f) ≡ φ(f). Since the friction force is not injective ζ is defined only in the part of the aging regime where the slope of the friction force is positive (cf. Fig. 1). Equation 5b shows that for ∂fζ/v̄ < 1 (frictionforce slope greater than unity) the diffusion coefficient becomes negative: D < 0. In this regime the system is driven towards the friction-force maximum via uphill diffusion. The 5 (September 18, 2000) requirement that the friction-force slope be greater than unity implies that the short-term response of the system be smaller than the steady-state response. Whether this condition is physically realizable (even in a system with memory) on a macroscopic scale is questionable [15]. Moreover, as the friction-force maximum is approachedD becomes singular (∂fζ → ∞). However, close to this maximum the adiabatic approximation is expected to break down as the time scale associated to diffusive relaxation of the control parameter vanishes. It is apparent that this diffusion-like instability, which is not present in the original BK model, contains many of the features previously attributed to SOC. It is characterized by an uphill diffusion and, in the adiabatic approximation, the diffusion constant becomes singular at the friction-force maximum. Either criterion has been identified as characteristic of SOC in continuum systems [9–11]. The nature of instabilities in Eqs. 3 is elucidated by performing a linear stability analysis of the uniform, steady-state solution e0 = 1− φ(v̄) , f0 = φ(v̄) , (6) to small perturbations [δe(0), δf(0)] exp(ωt+ikx). The roots of the characteristic polynomial are ω± = 1 ǫ2 ( μ± {μ − ǫ[1 + (1− φ 0 )k]} )